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We developed a new lattice Boltzmann method that al- 
lows the simulation of two-phase flow of viscoelastic liquid 
mixtures. We used this new method to simulate a bubble 
rising in a viscoelastic fluid and were able to reproduce the 
experimentally observed cusp at the trailing end of the bub- 
ble. 



I. INTRODUCTION 

The study of viscoelastic fluids is of great scientific in- 
terest and industrial relevance. Viscoelastic fluids are 
fluids that show not only a viscous flow response to an 
imposed stress, as do Newtonian fluids, but also an elas- 
tic response. Viscoelastic effects are almost universally 
observed in polymeric liquidstl, where they often domi- 
nate the flow behavior. They can also be obseijyed in 
simple fluids, especially in high frequency testing or in 
under-cooled liquidsB. Because most research into vis- 
coelastic liquids, especially that with an eye toward en- 
gineering applications, is directed toward polymeric liq- 
uids, the viscoelastic behavior of simple liquids is not as 
well known among researchers. The fact that the mani- 
festation of viscoelasticity does not require the presence 
of polymer molecules is at the heart of our approach, as 
will become clear in the description of the viscoelastic 
model. 

Although in most practical problems involving poly- 
meric materials the viscosities of the materials involved 
are so large that the creeping flow approximation is valid, 
the non-linearity introduced by the viscoelastic response 
of the liquid makes it difficult to treat any but the most 
simple cases analytically. In engineering applications the 
situation is often further complicated by the fact that 
the system is comprised of several immiscible or partially 
miscible components with different viscoelastic proper- 
ties. Examples of this include polymer blending, where 
two immiscible polymers are melted and mixed in an ex- 
truder, and the recovery of an oil- and- water mixture from 



porous bed rock. Simulation of these systems is very im- 
portant, but due to the complexities only few numerical 
approaches exist to date. Boundary element methods 
have been used to simulate such systems with varying 
degrees of success, but the allowable complexity of the 
interface morphology is very limited in such approaches. 
Lattice Boltzmann simulations have been shown to be 
very successful for Newtonian two-component systems 
with complex interfacedj, but for viscoelastic fluids 
lattice Boltzmann models, derived by Giraud et al.E 
are limited to one-component systems. 




FIG. 1. (a) An air bubble rising in Ivory® soap and (b) the 
simulation results for a low viscosity drop rising in a Oldroyd 
B fluid on a 256x1024 lattice for a drop of radius Rq — 35. The 
simulated bubble has a cusp that is rounded at the tip due to 
the finite thickness of the interface 3 lattice spacings). 

In this article we report the successful combination 
of both two-component and viscoelastic features into a 
two-dimensional lattice Boltzmann model. We used this 
model to simulate a bubble rising in a viscoelastic liquid 
(see Figure 1) and in this letter report the ffist successful 
simulation of the experimentally observed cusp. 



II. LATTICE BOLTZMANN 



* Centre for Systems Engineering and Applied Mechanics 
(CESAME), Universite catholique de Louvain, Batiment Eu- 
ler, Av. Georges Lemaitre 4, B-1348 Louvain- la- Neuve, 
Belgium 



We use a two-dimensional lattice Boltzmann model on 
a square lattice with a velocity set of {v^} = {(0,0), 
(0,1), (1,0), (0,-1), (-1,0), (1,1), (-1,1), (-1,-1), 
(1,-1)} and a corresponding set of densities {/i}, but 
following Giraud et alu we introduce two densities for 
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each non-zero velocity. We use a BGK lattice Boltzmann 
equation that contains the full collision matrix Ay 

fi{x + ViAt,t + At) = fi{x,t) 

+AtAij{f°{x,t)-fj{x,t)) (1) 

where the summation rule for repeated indices is imphed 
and the required properties of the equihbrium distribu- 
tions are discussed below. The local density is given 
by P = fi and the momentum by pu = fiVi. 

In order to simulate a two-component mixture we de- 
fine a second set of nine densities, {gi}, with an appropri- 
ate equilibrium distribution, {g^}. These densities rep- 
resent the density difference of the two components A 
and B as = X] = pA — pB^ where the total density 
introduced earlier is p = pA -\- pB- For the ^^s we choose 
a single relaxation time lattice Boltzmann equation 

gi {x + ViAt, t + At) = gi{x, t) 

+ — (2) 
r 

where r is the relaxation time and g^ is the equilibrium 
distribution. 

To use the lattice Boltzmann method in order to sim- 
ulate fiuid fiow, mass and momentum conservation have 
to be imposed. Mass and momentum conservation are 
equivalent to constraints on the equilibrium distribu- 
tions: 

i i i 

There will be further constraints on the permissible 
equilibrium distributions in order for the corresponding 
macroscopic equations to be isotropic and to simulate the 
systems in which we are interested. In the next two sub- 
sections we will summarize the physics that we want to 
incorporate and then we will discuss how it imposes con- 
straints on the equilibrium distributions and eigenvalues. 

A. Binary mixtures 

To simulate aJDinary mixture we follow the approach of 
Orlandini et al.u and begin with a free energy functional 
^ that consists of the free energy for two ideal gases and 
an interaction term as well as a non-local interface term: 

^[Pa,Pb]= / [TpaHpa)^TpbHpb) 

J X 

^>^PAPB + l^\d^{pA - pB)f] <^X, (4) 

where the densities pA and ps are functions of x. The 
repulsion of the two components is introduced in the A 
term and is a measure of the energetic penalty for an 
interface. When we write this free energy functional in 
terms of the total density, p, and the density difference, ^, 



we can derive the chemical potential, /i, and the pressure 
tensor, Pa/3, asB: 

/i = — = d^"^ - nd^d^cj), (5) 

^n{dcy,(t)df3(t) - -d^(l)d^(l)Sa(3 - (l)d^d^(l)Sa(3), (6) 

where 5 indicates a functional derivative and S^p is the 
Kronecker delta. For a two-component model wa^fix the 
further moments of the equilibrium distributiongj: 

E5N^ = '^"' (7) 

fiViaVifi = Pa/3 + pUaUff, (8) 

i 

^giViaVif3 = pSaf3 + (9) 
i 

Thus far, the model allows us to simulate a binary mix- 
ture that phase separates below a critical temperature of 
Tc = A/ 2. The surface tension, a, can be calculated an- 
alytically for a fiat equilibrium interface (j){y) orthogonal 
to the y direction as <j = z^: J^^{dy(l){y))'^ dy where the 
equilibrium density profile of <p also depends on k,. 



III. VISCOELASTICITY AND THE BOLTZMANN 
EQUATION 

Viscoelasticity was first proposed by Maxwell in his 
dynamic theory of gasesl. He used the simple argument 
that in the limit where there are no intermolecular colli- 
sions the fiuid in a container should behave like a solid: 

. . Then it can easily be shown that the pressures on the 
sides of the vessel due to the impacts of the molecules are 
perfectly independent of each other, so that the mass of 
moving molecules will behave, not like a fiuid, but like a 
solid." He goes on to deduce that the observed viscous 
behavior of fiuids is due to binary collisions that random- 
ize the directions of stress in the fiuid. Since the collisions 
are fast, but not instantaneous, the elastic properties of 
the fiuid are not completely lost, leading to the Maxwell 
model of viscoelasticity. 

Subsequently, derivations of hydrodynamics from the 
dynamic theory of gases have made the approximation 
of a purely viscous behavior because of the difficulties of 
deriving a continuum approach at the length scales of a 
mean free path of a molecule. In gases where lengths less 
than the mean free path are important kinetic theory for 
rarefied gases is used. 

There has recently been much activity in the research 
of the experimentally observed viscoelastic behavior of 
simple liquids that are undercooled. In this case, how- 
ever, viscoelasticity is not obtained because the relevant 
length scales were of the order of a mean free path, but 
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rather because of the correlations of subseauent cohisions 
as described in the mode couphng theoryB. 

The arguments of Maxweh, however, are stih vahd for 
describing the behavior of the Boltzmann equation, and 
viscoelastic properties can be derived from the Boltz- 
mann equation if the decay of viscous stresses is slow. 

The approach by Giraud et al. aims not at deriving 
a convected Maxwell fluid, but a convected Jeffreys fluid 
which is a mixture of a Maxwell fluid with a Newtonian 
fluid. A double set of densities is introduced allowing two 
stresses, one of which is chosen to relax quickly and is, 
therefore, a viscous stress, and the other, which is chosen 
to decay very slowly, represents a viscoelastic stress. The 
resulting model is a convected Jeffreys model that is often 
used to describe a polymeric fluid in a solvent. Care has 
to be taken for the choice of the collision matrix and 
the equilibrium distribution to ensure an isotropic model. 
The details of this one- component model are described in 
the publication by Giraud et aln 

A Chapman-Enskog expansion of the lattice Boltz- 
mann equations (|l|) and (g) gives the macroscopic equa- 
tions that our system simulates. Mass conservation gives 
the continuity equation: 

+ =0. (10) 

Momentum conservation gives a Navier Stokes equation: 

where the viscous stress is given by 

^af3 = ^oodpipUa) + ^oo 9^ {pU^)S a f3 - (12) 

The viscoelastic stress has the constitutive relation 

CTafB + ^Cr(l)a/3 = "(^0 - ^^oo) (9a(p^/3) + dfdipUa)) (13) 

where cr(i) represents the upper convected derivative of a. 
These equations are equivalent to the Navier- Stokes and 
Jeffreys equations only in the incompressible limit where 
da{pU(3) = pdaUf3. The fully compressible equations can 
only be simulated when a larger set of velocities is usedty. 
The conservation of the density difference leads to the 
convection diffusion equation 

(14) 

where is a diffusion constant given by = (r — 1/2). 

IV. SIMULATION OF A BUBBLE IN A 
VISCOELASTIC LIQUID 

We applied our method to a system similar to the ex- 
perimentally well-studied system of an air bubble rising 
in a viscoelastic fluid. In our simulation we represent the 



bubble using a phase-separated Newtonian drop of low 
viscosity in matrix which is viscoelastic by letting the re- 
laxation time in equation ( p!3|) depend smoothly on the 
density difference (j) between = 0.05 in the drop and 
= 66 in the surrounding fluid. We choose ^oo = 0.06, 
Uoo = 0.01, and uq = 0.01 in the drop and = 0.175 
outside. For the thermodynamic parameters we select 
T = 0.5, A = 1.1, and = 0.007, which corresponds to 
a surface tension of <j = 0.02. All units are in terms of 
the lattice spacings and the time steps At = Ax = 1. 
We introduce a forcing dependent on (f) so that the bub- 
ble is forced upward while the surrounding fluid is forced 
downward. We choose the total change in the momen- 
tum due to the forcing to be zero so that no walls are 
required in the simulations. 

We start the simulations without forcing and then pe- 
riodically increase the forcing after 10, 000 to 40, 000 iter- 
ations. We observe the change in the velocity u and store 
the distribution of (j) so that we have a way of judging the 
deformation of the bubble. 




(a) (b) 




(c) (d) 

FIG. 2. Shape for the simulated bubbles for different 
forcings in (lattice spacings )/(time step)^ (a) 1.6 • 10"^ (b) 
3.6 • 10"^ (c) 9.6 • 10"^ and (d) 1.96 • 10"^. In (c) a fit 
to the predicted form of the cusp (|x|^/^ is also shown. All 
simulations are after a 20,000 iterations at this forcing. 

Figure |^ shows the form of the drop for different forc- 
ings. At low forcings the drop is elongated in the flow 
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direction. This is in direct contrast to a bubble in a New- 
tonian fluid, which is flattened in the flow direction. At 
a larger forcing the bubble forms a cusp at the lower tip 
of the drop. For even larger forcings the drop starts to 
flatten in the flow direction. This sequence is in agree- 
ment with the experimental findingglil. Tha elongation 
of a rising bubble has been simulated befordlS, but this is 
the first time that the formation of a cusp has been sim- 
ulated. In Figure ^(c) it is shown that the cusp can be 
fitted-tp the functional form |x|2/3 predicted by Joseph 
et alx3 for a two-dimensional cusp created by the flow 
induced by two couter-rotating cylinders. 

Experimentally the formation of a cusp has been ob- 
served to coincide with a jump of nearly anjOpder of mag- 
nitude in the terminal velocity of the bubbletlO, although 
the mechanism remains disputed. On the one hand Bird 
et al. argue that surface-active impurities tend to immo- 
bilize bubble surfaces and hence retard the motion of gas 
bubbles. This discontinuous change in bubble shape may 
be responsible for the removal of the immirities, and thus 
lead to a jump in the final velocity. Liutll et al. alterna- 
tively suggest that the change in the shape of the bubble 
will make it more streamlined, and therefore increase the 
terminal velocity. 
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FIG. 3. Velocity, u, of the bubble in two different sets of 
simulations for a bubble of radius 35 in a 256x512 lattice. On 
the X-axis the iterations were multiplied by a scale factor so 
that corresponding forcings appear at the same point in the 
graph. No jump in the final velocity is observed at ~ 0.5 10^ 
iterations where the formation of a cusp is observed (indicated 
by arrow). The scale factors were 0.25 and 1. Velocity is 
measured in lattice units per time step. 

We examined the velocity for the rising drop as de- 
scribed above and found no jump of about half an order 
of magnitude as observed by Liu et al. in their exper- 
iment (see Figure |l|). Our simulations suggest that the 
jump in velocity they observe is not connected to a more 
streamline form of the bubble due to the cusp, but more 
likely to the presence of surfactants that are absent in 
our simulations. 



V. CONCLUSION 

We introduced a lattice Boltzmann model that can 
simulate viscoelastic two-component flows. We gave an 
intuitive explanation of the origin of viscoelasticity in 
our model and the model by Giraud et al. in terms of 
the original theory of MaxwelB. Simulations using this 
method have succeeded in reproducing the cusp at the 
end of a bubble rising in a viscoelastic medium that have 
eluded earlier numerical attempts with a more traditional 
boundary integral approach. 

The model has been successful in the qualitative sim- 
ulation of the bubble problem in two dimensions. We 
intend to extend the model to three dimensions in the 
future. This will also enable us to compare the results 
quantitatively with experiment. 
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